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Abstract 

The inelastic Lorentz gas in cooling states is studied. It is found that the 
inelastic Lorentz gas is localized and that the mean square displacement of 
the inelastic Lorentz gas obeys a power of a logarithmic function of time. It is 
also found that the scaled position distribution of the inelastic Lorentz gas has 
an exponential tail, while the distribution is close to the Gaussian near the 
peak. Using a random walk model, we derive an analytical expression of the 
mean square displacement as a function of time and the restitution coefficient, 
which well agrees with the data of our simulation. The exponential tail of the 
scaled position distribution function is also obtained by the method of steepest 
descent. 
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I. INTRODUCTION 



The statistical properties of dissipative systems have created considerable attention from 
the view point of non-equilibrium statistical mechanics. The analysis of gas kinetics has 
played an important role in the construction of statistical mechanics. Similarly, granular 
gases which are collections of inelastic particles are expected to be important to un- 
derstand non-equilibrium statistical mechanics of dissipative systems. To understand the 
statistical properties of dissipative systems, several attempts have been carried out on a 
collection of inelastic particles. 

The simplest model of granular gases is an inelastic hard core model. To analyze the 
time evolution of the velocity distribution function (VDF), we usually adopt the inelastic 
Boltzmann equation which is a generalization of the Boltzmann equation to the case of 
inelastic gases. There have been several investigations on VDF of the inelastic gas by the 
inelastic Boltzmann equation. |^-|6|j3|j7|-[T0| The approximate analytical solution of the in- 
elastic Boltzmann equation has been obtained by the method of expansion by the Sonine 
polynomial which is known as a fast convergent one and becomes standard technique to de- 
scribe both non-uniform elastic gases and uniform inelastic gas systems |^-|6|,|3|J^-[10| . They 
also suggested that the scaled VDF of the inelastic Boltzmann equation has an exponential 
high energy tail which has been confirmed by the direct simulation Monte Carlo 

method (DSMC) [§] or in molecular dynamics simulations of smooth hard inelastic disks 
in the homogeneous cooling state (HCS) [|J. On the other hand, some investigations on 
the spatial structure of the inelastic Boltzmann gas system have been performed [p JIT, 12 



Esipov and Poschel || have suggested that the diffusion length of particles obeys a power 



of a logarithmic function of time. Brey et al. |TT[ have derived the self-diffusion coefficient 
by means of a Chapmann-Enskog expansion in the density gradient of the particles in the 
inelastic Boltzmann gas system. Chen et al. [|1^] have shown in the molecular dynamics sim- 
ulation of sticky inelastic particles that the probability density function of inelastic particles 
obeys an exponential distribution at late stages. 



2 



There is, however, a strong restriction in the applicability of the inelastic Boltzmann 
equation and the method of the Sonine expansion. Inelastic collisions generate string-like 
clusters in a cooling state which sometimes reaches an inelastic collapse state. Thus, HCS 



cannot be kept in general []T3^ [T5| , [12[] . Thus, both of the inelastic Boltzmann equation and 
the method of the Sonine expansion which are only valid in dilute HCS [}3|,[5|,ffl[7|-[rH become 
no longer valid as time goes on |16|. 



On the other hand, the inelastic scatterings of classical electrons by ions can be modeled 
by the inelastic Lorentz gas fllTjillll in which few fraction of particles can move freely among 



many immobile scatters known as Drude's formula. Since the fixed scatters keep their 
initial configuration, the system can keep the state of dilute mobile particles if the initial 
configuration is dilute. Furthermore, there are some examples of the inelastic Lorentz gas 



in the natural phenomena such as the spontaneous interparticle percolation. [ll7| , |T8f The 



Lorentz gas is also familiar in dairy life such as the pedestrian motion among a cloud 



of persons and obstacles, the pachinko game and the pinball game |20| etc. There is an 
additional advantage that the analysis of the inelastic Lorentz gas model is easier than 
that of the inelastic Boltzmann equation because the inelastic Lorentz gas model treats the 
one-body problem: The situation of Lorentz gas is dilute mobile particles among almost 
immobile scatters, while Boltzmann g cLS IS cL collection of mobile particles. 

We believe that the systematic studies of the inelastic Lorentz gas system is also needed. 
The steady VDF of inelastic Lorentz gas and its percolation velocity under the uniform 



external field have been analyzed. |17|,[18| However, there have been no efforts to investigate 



the spatial structure of the inelastic Lorentz gas system without external field. It is important 
to discuss the spatial structure of the inelastic Lorentz gas without any external field as the 
first step of the study of the inelastic Lorentz gas system. Thus, we intend to investigate 
the spatial structure of the inelastic Lorentz gas in cooling state in this paper. 

The organization of this paper is as follows. In § 2, we introduce the model of our simu- 
lation. In § 3, we validate our simulation code by evaluating the mean square displacement, 
the velocity distribution and the position distribution in elastic collisions. We also explain 



the main results of our simulation and analysis for the mean square displacement and the 
position distribution of the inelastic Lorentz gas in § 3. We conclude our results in § 4. 

II. THE MODEL OF OUR SIMULATION 

In our simulation we adopt the event driven method for hard core collision of particles. 
2l|^3l The system consists of 40000 fixed smooth hard disks, that is 40000 scatters, and 



one tracer disk. The collision between the tracer and the scatters is inelastic hard core 
collision with a constant coefficient of restitution e. Through our simulation, we use dimen- 
sionless quantities in which the scales of length and velocity are respectively measured by 
the diameter of the disks d and the standard deviation of the initial velocity distribution 
U. Namely, both of all the scatters and the tracer are circular disks with their diameter 
d. The number density of the scatters in this system is 0.04d~ 2 , where the system size is 
L x = L y = lOOOd Because the area fraction is approximately 0.0314, the system of our 
simulation is dilute enough. We adopt the periodic boundary condition in this system. The 
system size is enough large to remove the boundary effect of this system, because the tracer 
is localized and have rarely gone across the boundaries. We average each result over 10000 
samples which consist of 100 different initial configurations of the scatters and 100 different 
initial velocities of the tracer. Each initial configuration of the scatters is at random. Each 
tracer starts from the center of the system and its distribution of the initial velocities is 
made from normal random numbers, where its average is zero and its standard deviation is 
U. Each sample has been calculated until 100 inelastic collisions. 

In each inelastic collision, the velocity of a tracer tangential to the collision plane are 
preserved, while the velocity normal to it changes with the restitution coefficient e which 
is assumed to be a constant for all the collisions. The precollision velocity of the tracer v* 
yields the postcollision velocity of it v with v • a — — e(v* • a), where a is a unit vector 
pointing from the center of the scatter to the center of the tracer ( Fig.|T] ) . Thus, the 
expression is given by 
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v = v*-(l + e)(v*-a)<T. (1) 

This relation between the precollision velocity and the postcollision velocity in the inelastic 
Lorentz gas system is quite different from that in the inelastic Boltzmann gas system. 

In the simulation of the inelastic Lorentz gas, the velocity of the tracer quickly decreases. 
For the sake of long time simulation, we adopt the rescaling of the speed of the tracer as 
follows. After each collision, the magnitude of the postcollision velocity |v| has been put to 
that of the precollision velocity |v*| , i.e. the magnitude of the initial velocity of the tracer 
has been preserved for all the collisions in each sample. Note that we use non-scaled values 
of v the magnitude of which decreases by collisions for our analysis. 

III. RESULTS: SIMULATION AND THEORY 

Now, we explain the main results of our simulation and analysis for inelastic Lorentz gas. 
We study the mean square displacement and the position distribution of the tracer. 

A. Check of our simulation code 

Our simulation code has been validated by the test simulation in which we have evaluated 
the mean square displacement, the velocity distribution and the position distribution of the 
tracer in the elastic collisions, i.e. e = 1. For elastic scatters, we expect that the tracer 
displays random walk motion when the initial speeds of the tracer are identical and that the 
Gaussian is the stable velocity distribution. 

Corresponding to a random walk model, we have fixed the magnitudes of 100 different 
initial velocities of the tracer, i.e. |v | = U. In each sample, we calculate the motion of tracer 
until 10000 elastic collisions. The other conditions of the test simulation are same as those 
of the simulation in inelastic collisions. We have checked that the mean square displacement 
of the tracer (r 2 (t)) is almost a linear function of time, where (• • •) indicates the average 
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over the 10000 samples. We also confirm that all the data of the position distribution g(r, t) 
can be scaled as 

g{r,t) = -=L=G{R), (2) 

r ¥W) 



where G(R) is the scaled position distribution function and R = r/^(r 2 (t)} is the scaled 
distance from the center of the system at which the tracer started. We check that the 
scaled position distribution G(R) is almost on the Gaussian curve as expected from the 
conventional random walk. Furthermore, the velocity distribution of the tracer becomes 
almost the Gaussian. 



B. Mean square displacement 

Let us discuss the statistical properties of inelastic Lorentz gas. At first, we evaluate 
the time evolution of the thermal velocity and discuss the mean square displacement in the 
inelastic system. 

In the inelastic Lorentz gas system, many inelastic collisions happen at random. The 
average number of collision times per unit time is evaluated as 2dvo(t)po, where d is the 
diameter of the tracer and the scatters and po is the number density of the scatters. The 
thermal velocity vo(t) is defined as 

v (t) = ^/vV(v,f)dv, (3) 

with the velocity distribution of the tracer f(v,t) and the average density of the tracer 
n = j f(v,t)dv. Therefore the temperature of the system is denoted as 

T(t) = \mvl{tl (4) 

where m is the mass of the tracer. ||,|,[|,[3] 

In each collision, the energy loss of the tracer is considered to be proportional to the 
temperature of the system T(t). Thus, we can estimate the energy loss of the tracer per 
unit time as 
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= -2dv (t)p h(e)T(t), (5) 

where h(e) is a dimensionless function of the restitution coefficient e. h(e) represents a 
cooling rate of the inelastic Lorentz gas system. Substitution of eq. (fD into eq. (^j leads to 

^ = -rfpoMe)^W- (6) 
Solving eq. (|]) , we obtain the thermal velocity t>o(t) as 

= — — , w u > (7) 



where t>o(0) is the initial thermal velocity. In the inelastic Boltzmann gas system, Haff [24 
has derived the same expression as that of eq.(^), where the cooling rate is assumed to be 
proportional to 1 — e 2 . 

Here, we can deduce the expression of h(e) from our simulation as follows. Figure |2| 
shows (vo\t) - v \0))/(l - e)U versus tUd 1 plots for several values of the restitution 
coefficient, namely, for e = 0.4, 0.6, 0.7 and 0.9. All the data are on the same straight line 
which is also shown in Fig.|2|. Corresponding to the results of our simulation, we set the 
number density of the scatters po — 0.04c? -2 and the initial thermal velocity t>o(0) = \/2U 
from eq.(|). Thus, Fig.|] indicates that h(e) is a linear function of the restitution coefficient 
e as 

h{e) ~ 1.45 x (1 -e), (8) 

whereas Martin and Piasecki |18| have assumed energy losses at collisions in inelastic Lorentz 
gas system are proportional to 1 — e 2 . This Haff's law, h(e) ~ 1 — e 2 , can be derived from 
the Gaussian VDF. We will demonstrate that the cooling rate of the inelastic Lorentz gas 
system is not proportional to 1 — e 2 and close to eq.(H) in terms of non-Gaussian nature of 
VDF as follows. 

The time evolution of VDF for inelastic Lorentz gas without any external field is in 
general described by a rate equation. For elastic gases, the rate equation is known as the 



Boltzmann equation. Similarly we can derive a rate equation of VDF for inelastic Lorentz 
gas systems (see Appendix A), i.e. the inelastic Lorentz gas equation : 



dt 



dp J (v^)[^-/(v,t)]da. (9) 
Jv-ct>o e 



We assume that eq.(0) admits an scaling solution depending on time only through the 
thermal velocity v (t) defined in eq.(|) as 

/(v,i) = -^-F(c), (10) 

with the scaled velocity of a tracer c = v/vo(t) and the average density of the tracer 
n = J f(v,t)dv. F(c) is the scaled VDF of the tracer which is independent of time. Here 
we define the dimensionless collision integral as 

1(F) =f (ca)[^-F(c)]da, (11) 

and its moments as 

M p = - J Vj(F)dc. (12) 

Multiplying eq.(Q) by mv 2 /2 and integrating by v, we obtain the rate of change of the 
temperature from eqs. @, (0), and ([12]) as 

dT(t) 



dt 

Comparing eq.([L3|) with eq. (0), we obtain 



dv {t)p M 2 T{t). (13) 



He) = ^. (14) 



Both sides of eq.(|T^) represent the cooling rate of the inelastic Lorentz gas system. When 
we adopt the lowest order Sonine expansion in eq.([B3|), the expression of M 2 as a function 
of the restitution coefficient e is given as 

65 15 1 

M 2 = 0F(1 - e 2 )(— + — e 2 - — V-383 + 828e 2 + 900e 4 ). (15) 
o4 62 o4 



The details of the derivation of eg. ( |B27|) are written in Appendix B. Figure |3| shows the both 
sides of eq. (fH|) , i.e. M 2 /2 from eq.( [B27| ) and h(e) from eq.([8]) which has been estimated 
from the data of our simulation, as a function of the restitution coefficient e. There are 
no significant differences between M2/2 and h(e), though M2/2 does not exist in the high 
inelastic region. This means that Haff's law, h(e) ~ 1 — e 2 , may be replaced by eq.(^J) in 
the inelastic Lorentz gas system as a result of non-Gaussian nature. 

Now, let us discuss the statistical properties of the inelastic Lorentz gas based on a 
random walk model. To simplify the expressions, we introduce a and (3 as, 

a = f3 = dp h(e). (16) 

Using a and /3, the average time interval tn from iV — 1th collision to iVth collision between 
the tracer and the scatters is given as, 

T N ~—t — - = l(a + pt N - 1 ), (17) 
«o(*jv-i) 

where / denotes the mean free path of the tracer, and tjv-i is the total time until the N — 1th 
collision. From the relation t^-i = t~\ + t 2 + • • • + tjv-i, we obtain as a function of the 
number of collision times N, 

TN = la{\ + l(3) N -\ (18) 

Therefore the total time until the iVth collision is given by 

t N = Tl +r 2 + • - • + r N = ^{(1 + l(3) N - 1}. (19) 

Rearranging eq. (p!9|) , the number of collision is obtained as a function of the total time 
until the iVth collision as 

lnf^tjv + 1) 

N = a -. (20) 

ln(l + 1(3) 1 J 

Note that eq.(^) holds even for the elastic limit 

N _, f JL for e -»• 1. (21) 



Here we have used the result of eq.(U) , i.e. h(e) — >• for e — ► 1. 

If we assume that the inelastic Lorentz gas system corresponds to the random walk 
system with the reduction of hopping velocity, the effect of inelasticities can be absorbed in 
the number of collisions N. Thus, we deduce the expression of mean square displacement of 
a tracer as 

^ = eN - ' 2 WW)' (22) 

where we use eq.(^0]) and the suffix e of the bracket emphasizes the quantity with the 
restitution coefficient e. Let us evaluate some physical quantities in eq. fl2~2|) . Using the 
values of our simulation, we obtain a = l/\/2U~ l , (5 ~ 0.058 x (1 — e)d~ l from eqs.(§) and 
(|16|). We set the mean free path of the tracer as / ~ 13.5<i instead of the theoretical mean 
free path I = 12.5<i. 

Substituting these physical quantities and eq.(||) into eq.(|22|), we obtain the mean square 
displacement as a function of time and restitution coefficient. Figure ^ shows {r 2 (t)) e x ln(l + 
lf3)d~ 2 versus (1 — (^iUdr 1 plots for e = 0.4, 0.6, 0.7 and 0.9. All the results of our simulation 
shown as the symbols are on the same curve. Our results clearly demonstrate that the tracer 
is localized in the case e < 1 and that the mean square displacement seems to obey a power 



of the logarithmic function of time. It should be emphasized that eq.(]22|), shown in Fig.|], 
becomes identical with the results of our simulation as time goes on. The reason why eq.(p2|) 
deviates a little from the data of our simulation in the early time region in Fig.^ is that few 
collisions have happened in this region. 

C. Position distribution 

In the inelastic Lorentz gas system, the position distribution function of the tracer g e (r, t) 
for the restitution coefficient e is scaled by the mean square displacement (r 2 (t)) e . The scaled 
position distribution function of the tracer G(R) is defined as 



G{R) = \] (r 2 {t)) e x g e (r,t), (23) 
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where R = r/y(r 2 (t)) e is the scaled distance from the center of the system at which 
the tracer started. G(R) is independent of the restitution coefficient e and time. Let us 
study the scaled position distribution of the tracer G(R). Figure [5] shows the results of 
our simulation and analysis. The symbols are results of our simulation which show the 
scaled position distribution of the tracer for e = 0.4, 0.6, 0.7 and 0.9 at four different times 
t = 3000dU-\ 300dU-\ bOOdU" 1 and lOOOdU' 1 , respectively. All the data of the scaled po- 
sition distribution G(R) are almost on the universal curve. It should be noted that the curve 
is close to a Gaussian near its peak, while the tail seems to obey an exponential function. 



This exponential tail has been observed in a simulation of granular gas system. [[12 

On the other hand, the broken line and the dash-dotted line in Fig.|5| are results of our 
theoretical analysis which will be explained later. The broken line corresponds to the data 
in the tail region of the scaled position distribution, while the dash-dotted line corresponds 
to the data in the peak region of the scaled position distribution. Theoretical curves well 
explain the result of our simulation. Now, let us explain how to obtain the theoretical curves. 

At first, we assume that the tracer obeys a renormalized diffusive motion as assumed 
in the previous subsection. This assumption is validated by the derivation of the diffusion 
equation of the tracer by means of a Chapman-Enskog expansion in the density gradient 
of the tracer. [|i~0[|l"ll The time evolution of the position-velocity distribution of the tracer 
p(r, v, t) is given by the inelastic Lorentz gas equation (see Appendix A): 

M^M + v . Vj 9(r, v, t) = Ap(r, v, t), (24) 

where the operator A is defined as 

Ap(r, v, t) = d Po [ (v • a) [ P ^ V *^ - p ( r , v , t)]da. (25) 
Jv-s->o e z 



Integrating this equation by v, we obtain 

dn(r, t) 



Of 

where the density of the tracer n(r, t) is defined as 
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+ V-J(r,t) = 0, (26) 



n(r,t) = J dvp(r,v,t), (27) 
and the flux of the tracer J(r, t) is defined as 

J(r,t) = /dvvp(r,v, i). (28) 



Note that the integral of the right hand side of eq. (|24[) by the velocity of the tracer becomes 
zero, because it denotes the effect of change in the velocity of the tracer. 

Here we adopt the Chapman-Enskog method in which the position-velocity distribution 
of the tracer p(r, v, t) is written as a series expansion in a formal parameter e, 

p(r, v, t) = p {0) (r, v, t) + ep {1) (r, v, t) + e 2 p {2) (r, v, t) + • • • . (29) 



e implies an implicit Vn(r,t) which provokes diffusion of the tracer. Consistently, eq.(|26|) is 
expanded as 

^M = £ e ^dPn(v,t), (30) 



3=0 



with 



d[ j) n{v,t) = -V (0) • J (i) (r,t), (31) 



and 



jW(r,t) = y dwp^fov,*). (32) 



Note that 



V = e V (0) + e 2 V (1) + ■ ■ • . (33) 
The zeroth order of eq. (f2~4"D , thus, becomes 



OTCQ flp(°)(r,v,t) = (0) 



(34) 



A comparison of this with eq.(^j) indicates that p^(r, v, t) is proportional to /(v, t). Owing 
to the normalization, 
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P (0) (r,v,t) = ^/(v,t). (35) 



?? 



Using this relation and eq. (|32"D , j(°/(r, t) = 0, i.e. the flux of the tracer vanishes to zeroth 
order, which leads to d^n(r,t) = 0. Additionally, from eq.(|35|), dp(°\r,v,t)/dn(r,t) = 
f(v,t)/n. To the first order in the density gradient e, eq.(|24]) leads to 



{dt n{r ' t)} dn(r,t) WW~ 

+ (v • V (0) n(r,t))^^ = ApW(r, v,t). (36) 
Using the results of the zeroth order and eq.([5D, this equation becomes 

mT{t) &k) +A)p(1)(r ' V)t) = (V ' V( °^ (r,t)) /(v,*), (37) 
with = 2dvo(t) poh(e) . The solution p^(r, v, t) is proportional to V^n(r,t) : 

p( 1 )(r,v,t) = B(v,T(t))-V(°)n(r,t), (38) 

where the vector B(v, T(t)) also depends on t through T(t). Substitution of this into eq.(|37D 
leads to 

mT{t) Wif) + A)B(v ' T{t)) = ^ /(v ' t] - (39) 

Here, the contribution to the particle flux of first order in the density gradient is given by 

J«(r,t) = - J D e (t)V (0) n(r,t), (40) 

with 

D e (t) = ~ J dvv-B(v, T(t)). (41) 

Here we have used symmetry considerations in two- dimension. Using the diffusion coefficient 
D e (t), eq.(4l]) becomes 



mnt)-^--u D )D e (t) = --, (42) 

ol (t) m 
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with 



/(AB)-vdv 

"° S ~ JB- vdv • (43) 



Because £(t) ~ T^if), dimensional analysis requires D e (t) ~ T*(t). Therefore, 

3£> e (t) _ geg) 

#r(t) 2T(t) ' 

Substitution of this into eg. ([42]) becomes 



(44) 



D , (t) = m (VB _ m r , (45) 

From eqs. fliOD and (|3lD , we obtain the diffusion equation of the tracer, 

dPn(r, t) = D e (t)V {0)2 n(r, t), (46) 

with the diffusion coefficient D e (t) in eg. fl45|) . The expression of the diffusion coefficient 
D e (t) is obtained explicitly in eg. (|C7|) with the aid of eg. ( |B18| ) by the method of the Sonine 
expansion, though the diffusion coefficient D e {t) in eg. (P7|) gives a little deviated results 
from the results of our simulation in the inelastic Lorentz gas system as shown in Appendix 
C. Now the random walk picture is validated by this standard method. 

Therefore, the tracer obeys a renormalized diffusive motion; the original position distri- 
bution function of the tracer g e (r,t) is expected to be given by 

1 r 2 
^Io D e{y,t')dt' 8 J D e (V, t')dt' 

Here we have taken into account the distribution of hopping velocity; we have adopted the 
diffusion coefficient as 

(48) 

where k is undetermined dimensionless coefficient which will not be indicated explicitly, r 
is time interval between two adjacent collisions, 1 is the hopping vector of the tracer, and V 
is the magnitude of the scaled hopping velocity defined as 
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V 



\[r_ 

v (t') 



(49) 



The expression of eq.(^HJ) is derived in Appendix D. Introducing 77 as 



71 V*{r*(t))e ' (50) 



we obtain the expression of r\ from eq. (|22|) as 



4/dn(l + Z/3) 

r/ = ^ , (51) 



which is independent of time. Thus, from eqs. fl47|) and (pOj) , we can derive the original scaled 
position distribution function of the tracer G(R, V) as 

G(R,V) = -^Re'S* . (52) 

Here we assume the scaled hopping velocity corresponds to the scaled velocity of the 
tracer c = v/t> (t). In our simulation, all the data of the scaled velocity distribution of the 
tracer F(c) can be fitted by the Maxwell-Boltzmann distribution 7r _1 e _c2 as Fig.[| shows, 
though the high energy tail of the tracer deviates from the Gaussian (see Appendix B). 
Therefore, we assume the distribution of the scaled hopping velocity is almost Gaussian. 
Now, the scaled position distribution function of the tracer G(R) is expressed by 

roo p— V 2 1 n 2 9 poo 

G{R)= —Re~^2nVdV = -R / e~ p(v ' R)R dV. (53) 

JO IX T] V T] JO 

The second equality in eq . (|53"D defines p(V, R) as 

\nV V 2 1 
P{V ,R) = — + - + —. (54) 

We analyze eq.(p3|) by the method of steepest descent in order to obtain the tail region 
of the scaled position distribution function. Designating V which satisfies dp(V, R)/dV = 
as V*(R), the expression of V*(R) is 



V*(R) 



1 V^TSW 

\ 4 + 4^ • (55) 
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By the method of steepest descent, from eq.(^) we can derive the approximate G(R), 

i.e. 



2 r°° 
G(R) ~ -R / exp 
n Jo 



2 '■ x 
77 Jo 

2 / -p(V*,R)R 2 



-p(v*,r)r 2 - v 2 7 ^; q y=y ,fl 



'/ V ^ 9V 2 |vw* 
Equation fl56f ) is reduced to 

G(R) 



dV 

e-PW'W. (56) 



7ri? -Jzr 



c 



(57) 



in the limit of R — > 00. Thus, we have shown that the tail of G(i?) obeys an exponential 
law. Now we use 

g(*> = 2 i Lsw\ (58) 

' V 1 dV 2 \V=V* 

instead of eg . (|56|) . Here 7 is a fitting parameter. We estimate 77 ~ 0.25 and 7 ~ 7.39 from 
the comparison of eq . (|58|) and the data of our simulation on the scaled position distribution 
in the tail region. 

In order to obtain the peak region of the scaled position distribution function, we replace 

V in eq.fl5"2|) by constant V. Thus, eq.(^) becomes 

G(R) ~ ~Re~S*. (59) 
r]V z 

Thus, we expect the Gaussian near the peak. Adopting the most probable V as V, i.e. 

V = l/v2j V is estimated to be 0.75. 

In Fig|| eq.(^) with r\ = 0.25 and 7 = 7.39 is expressed by the broken line, while eq.(|59"|) 
with rj = 0.75 and V = l/y/2 is expressed by the dash-dotted line. Figure [| shows that 



eq. (p8|) is in good agreement with the simulation data of the scaled position distribution in 
the tail region, while the data of our simulation on the scaled position distribution in the 
peak region are well-fitted by eq. (|59"D . 
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IV. DISCUSSION 



Some approximate analytical solutions of the inelastic Boltzmann equation obtained by 
the method of the Sonine expansion, namely scaled VDF of a granular gas, are in good 
agreement with the results of DSMC in which HCS holds [|3],0|| or of the molecular dy- 
namics simulations of smooth hard inelastic disks in HCS [Q. However, it has been also 
demonstrated that the Sonine expansion diverges if the inelasticity is high in the inelastic 
Boltzmann gas system. B Solutions by the Sonine expansion cannot express the definite 
exponential form in the high energy tail, because the expression of eq. (|B4j) does not have 
terms of odd power of c, while the definite exponential form contains odd power of c. This 
means that we cannot express the scaled velocity distribution of the inelastic Boltzmann 
gas system of high inelasticity by eq . (|B4j ) even if we obtain the perfect form of eq.(|B4]). 
The reason why a 2 does not exist in the inelastic Lorentz gas system of high inelasticity (cf. 
eq. (p26|) ) is thought to be the same reason mentioned above: the solution by the Sonine 
expansion does not have terms of odd power of c and it cannot express the scaled velocity 
distribution of the inelastic Lorentz gas system of high inelasticity because its high energy 
tail is definitely exponential. Furthermore, in both of the inelastic Boltzmann gas system 
and the inelastic Lorentz gas system, some of the scaled velocity distributions obtained by 
the method of the Sonine expansion become negative for large c, which is obviously an odd 
result. 

We assume the scaled VDF of the tracer obeys the Gaussian distribution in § 3.3, while 
we have derived the non-Gaussian scaled VDF of the tracer in Appendix B. Though Fig. 
6 shows the scaled VDF of the tracer seems to obey the Gaussian, its high energy tail is 
thought to become exponential. The Gaussian scaled VDF shown in Fig. 6 is rationalized by 
the insufficiency of the number of the samples in our simulation. The non-Gaussianity of the 
scaled VDF is not essential for the exponential tail of the scaled position distribution function 
of the tracer, while it is essential for the fact that the cooling rate is not proportional to 1 — e 2 
but close to 1 — e. However, the cooling rate is almost proportional to 1 — e 2 in the inelastic 
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Boltzmann system in which the scaled VDF is non-Gaussian. The difference is considered to 
be due to the difference between one-body problem and many-body problem, i.e. the collision 
relation between precollision velocity and postcollision velocity of the inelastic Lorentz gas 
system is quite different from that of the inelastic Boltzmann gas system and the mean free 
time of the inelastic Lorentz gas system also shows different dependence on time from that 
of the inelastic Boltzmann gas system. 

We have analyzed the spatial structure of the inelastic Lorentz gas system by both of the 
method of the Sonine expansion and the random walk model. The mean square displacement 
of the tracer by the random walk model is identical with the results of our simulation, while 
the mean square displacement by the method of the Sonine expansion does not agree well 
with them. Furthermore, the exponential tail of the scaled position distribution function of 
the tracer cannot be derived by the method of the Sonine expansion. Thus, to analyze the 
inelastic Lorentz gas system, the random walk model is considered to work better than the 
method of the Sonine expansion, though the random walk model also has some problems, 
e.g. we have used the mean free path of the tracer which is a little different from the 
theoretical mean free path. 

V. CONCLUSIONS 

We have found that the mean square displacement of the inelastic Lorentz gas obeys 
a power of the logarithmic function of time and that the scaled position distribution of 
the inelastic Lorentz gas has an exponential tail, while the peak is the Gaussian. We have 
analytically derived the mean square displacement as a function of time and restitution 
coefficient which is identical with the data of our simulation. We also analytically obtained 
the scaled position distribution function the tail of which is exponential and the peak of 
which is the Gaussian. 
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APPENDIX A: TWO-DIMENSIONAL INELASTIC LORENTZ GAS EQUATION 

We explain the derivation of two-dimensional inelastic Lorentz gas equation which can 
not be derived directly from the inelastic Boltzmann equation. 

We designate the position-velocity distribution function of the tracer at time t as 
p(r, v, t). The change of p(r,v, t) during minute time 5t is expressed as 

[p(r, v, t + 5t)- p(r, v, t)]drdv = -r_ + T + + T floW} (Al) 

where r + is the gain term of p(r, v, t) by collisions during the minute time, while T_ is the 
loss term of p(r, v, t) by collisions during the minute time. Tfl ow indicates the change in 
p(r, v, t) due to the flow into or out of dr, free- streaming term, which is given by 

r flow = -drdv5t(v • V)p(r, v, t). (A2) 

T_ and T + are defined as 

r_ = J p( T , v,t)drdvp \v\Stdb, (A3) 

and 

r + = J p(r,v*,t)drdv*/9o|v*|<Jtd& = J p{r,v*,t)dr^-p \v*\6tdb, (A4) 

respectively. Here v* denotes precollision velocity of the tracer, p is the mean number 
density of scatters, e is the restitution coefficient, b indicates the impact parameter ( see 
Fig.|7| ). Note that we use the Jacobian dv* = dv/e. Substitution of eqs.(|X3]),(|X^) and ( |A2| ) 
into eq.flAj]) leads to 
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+ (v . V )p(r, v, t) = po J [[v*l P(r 'J V) - | v|p(r, v, t)]d&. (A5) 
It is often convenient to rewrite the integral over the impact parameter as the unit vector 



a as 



|v|d6 = B(\v\,a)da, (A6) 

Here da = d(7r — if>) and b = dsin(iT — ip) ( see Fig.[7| ). d is the diameter of the tracer and 
the scatters. Thus, the expressions of -B(|v|, 5") and B(\v*\,a) are given by 

B(\v\,a) = |v|dcos(7r -ip) = d(v ■ a), (A7) 

and 

B(\v*\,a) = |v*|rfcos(7r -ip) = -d(v* ■ a) = di^— ^, (A8) 



respectively ( see Fig.|7] ). From eqs.( |A5| ), (|A6|) , ( |A7| ) and (|A8|), we obtain 



^ (r ' V>t) + (v • V)p(r, v, t) = rfp / (v • a) [^1A _ p(r , v , f) ]d *. (A9) 
at Jv&>o e z 

This is nothing but the two-dimensional inelastic Lorentz gas equation. Integrating this 
equation by r, we obtain eq.@. 

APPENDIX B: COOLING RATE M 2 



We derive the expression of M 2 by the method of Sonine expansion. 
From eq.(f|), eq. ([l3|) is rearranged as 

dv (t) Po dM 2 v$(t) 



(Bl) 



dt 2 

Using eqs.(p|) ,(pT3|) and ( |B1|) , we obtain the scaled equation for the inelastic Lorentz gas as 



J(F) = ^(2 + c~)F(c). (B2) 



Integration of the product of c p and eq.( |B2|) by c makes 
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M P = ^p(cf), (B3) 

where (• • •) denotes an average over F(c). 

In the limit of small dissipation, F(c) approaches a Gaussian, i.e. F(c) ~ 0(c) = 
7T _1 exp(— c 2 ). Therefore, in the low inelasticity region, a systematic approximation of F(c) 
may be expressed by the Sonine expansion H,H,[7EfflII 



i.e. 



F(c)=0(c){l + ^aA(c 2 )}. (B4) 
P =i 



The first few two-dimensional Sonine polynomials are 

S (x) = 1 (B5) 
S 1 (x) = -x + l (B6) 
S 2 (x) = ~ x 2 -2x + l. (B7) 

It should be noted that Sonine polynomials satisfy the orthogonality relations 

j 'dc4 > (c)S p (c 2 )S p ,(c 2 )=S pp ,, (B8) 

where S pp i is the Kronecker delta. From these orthogonality relations, we have 

a p = J dc^(c 2 )F(c) = (S p (c 2 )). (B9) 

From eqs. (|B6|) , (p7|) and (p9|) , a x and a 2 is given as 

a 1 = -( c 2 ) + l = (BIO) 
a 2 = ^-2(c 2 ) + l = ^-l, (Bll) 

respectively. Here we have used the fact that substitution of p = 2 into eq.([B3|) leads to 
(c 2 ) = 1. Substitution of p = 4 into eq.(p3|) leading to 



M 4 = 2M 2 (c 4 ), (B12) 

we obtain 
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4M 2 (a 2 + l)-M 4 = 0. (B13) 
On the other hand, the moments of the dimensionless collision integral M p is calculated 



as 



M p 



[dc [ da(c ■ a)F(c)(c p - c' p ), (B14) 
J Jc-a>o 



where d indicates the postcollision velocity of the tracer. From eq.( [B13| ), we need the 



expression of M 2 and M 4 to obtain a 2 . In general, it is impossible to obtain such the moments 
without the complete information of F(c). Thus, we adopt an approximate method to obtain 
a closed relation for M 2 , M 4 and a 2 , i.e. we approximate the scaled velocity distribution 
function as F(c) = 0(c){l + a 2 S' 2 (c 2 )} using the expression of eq. flB4] ) with a p = for p > 3. 
Adopting such the approximation, the second moment is calculated as 

M 2 = v^(l-e 2 )(l + |a 2 ). (B15) 

o 

The forth moment is calculated similarly as 

M 4 = v^(l - e 2 )(3 + 2e 2 )(l + ^a 2 ). (B16) 

o 

Inserting these moments into eq.( pi3| ), we obtain 



a 2 = — (1 + 30e 2 ± V-383 + 828e 2 + 900e 4 ) if 1. (B17) 



We adopt the stable solution of a 2 0, i.e. 



1 



a 2 = — (1 + 30e 2 - V-383 + 828e 2 + 900e 4 ), (B18) 

for the inelastic case. Another solution of a 2 in eq.( |B17| ) is unstable. 

Here we discuss the stability of the coefficient a 2 in eq.( |B17| ). To analyze the stability of 



a 2 we write eq. ( ID ) in a more general form 



71 

f(v,t) = ^-F(c,t), (B19) 
v o{ 1 ) 



which leads to 
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dpov (t) dt 2 

in essentially the same way as eg. ([B"2"|) . Using F(c,t) = 0(c){l + a 2 (t)S 2 (c 2 )} and (c 2 ) = 1 
which is bound, substitution of p = 4 into eq.( |B20| ) leads to 

^ = fiP0|(i) {4M2(a2(i)+1) _ Ml}i (B21) 

for the coefficient a 2 (t). To discuss the linear stability of the coefficient a 2 (t), we divide it 
into the equilibrium part a 2 and the time-dependent part 5a 2 (t), i.e. 

a 2 {t) = a 2 + 5a 2 (t). (B22) 



Note that a 2 is a 2 in eq.( |Bf 8| ). Substituting eq. ( |B22|) into eq.( B21|) and ignoring the higher 
order 0(5a\{t)) terms, we obtain 

^ = d fV^(l - e*)H(e)v (t)6a 2 (t), (B23) 
which can be solved for a 2 (t) as 

.dp 



5a 2 (t) = 6a 2 (0) exp[-^V^(l " e 2 )H(e) f v (t)dt], (B24) 
2 Jo 



where 



H(e) = 3a 2 -( 1 ^+ 1 -). (B25) 

Thus, the stability of the coefficient a 2 {t) depends on the signs of H(e), i.e. a 2 (t) is stable 
if H (e) is negative. Figure [8] shows H(e) as a function of the restitution coefficient e. Here 
we substitute the expressions of a 2 in eq.( |B17|) as a 2 . From Fig. ||, it is shown that a 2 in 



eq. (pi8|) is stable, while another a 2 in eq.( [B17|) is unstable. 



It should be also mentioned that the elastic case e = 1 is a singular point in eq. (pi8|) , i.e. 
F(c) = </>(c){l + a 2 S 2 (c 2 )} does not become the Gaussian even if e = 1, and that eq. ( |B18| ) 
has real values only in the region 
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The reason why a 2 does not exist in the inelastic Lorentz gas system of high inelasticity 
is thought to be that the solution by the Sonine expansion does not have terms of odd 
power of c and it can not express the scaled velocity distribution of the inelastic Lorentz 
gas system of high inelasticity by eq. (|B4|) even if we obtain the perfect form of eq. ( |B4| ) 



because its high energy tail is definitely exponential. It has been also demonstrated that the 
Sonine expansion diverges if the inelasticity is high in the inelastic Boltzmann gas system. 
0] This means that we also can not express the scaled velocity distribution of the inelastic 
Boltzmann gas system of high inelasticity. 

From eqs.( |B15l) and ( B18|) , the expression of M 2 as a function of the restitution coefficient 



e is given as 



,65 15 



M 2 = ^(1 - e 2 )(— + —e 2 v / -383 + 828e 2 + 900e 4 ). (B27) 

v v A 64 32 64 ; v ; 



Judging from the expression of eq.( B15| ), the reason why M 2 is proportional not to 1 — e 



2 



but to 1 — e is due to the deviation of the scaled velocity distribution function from the 
Gaussian, i.e. if a 2 = in F(c) = </>(c){l + a 2 S 2 (c)}, M 2 in eq.( [B15|) becomes proportional 



to 1 — e 2 and the cooling rate of the inelastic Lorentz gas system becomes proportional to 
1 — e 2 analytically. Note that because we multiply a 2 by 1 — e 2 in eq. (pi5|) , M 2 holds even 
at the elastic case, i.e. e = 1, though eq. (|B18j ) does not hold at the elastic case. 



APPENDIX C: DIFFUSION COEFFICIENT D E (T) 

We derive the expression of the diffusion coefficient D e {t) in eq.(^) explicitly by the 
method of the Sonine expansion. 

Calculating eq.(f|) in a similar way as eq.( [B14|) , we obtain 



f dvF(v)ylX(v) = d Po [ dv f da(v • a)X(v)[Y(v') - Y(y)], (CI) 

where v' is the postcollision velocity of the tracer: 

v' = v- (l + e)(va)a. (C2) 
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Here we approximate B(v, T(t)) in the first Sonine approximation as 

B(v,T(t))~v/ a2 (v,t), (C3) 
where the velocity distribution function of the tracer f(v,t) is approximated as 

f(v,t) ~ f^t) (C4) 



exp[-— — v 2 ]{l + a 2 (l - — -v 2 + ——v 4 )}. 



2vrT(t) 2T(t) JL v T(t) 8T 2 (t) 



(C5) 



Use of eq. flClD into eq.(^) leads to 



dpo(l + e) gTrTg) , 3 

" D = 2 V ~^ (2 + I a2) - (C6) 



Substitution of this expression of vn into eq.(f4"5|) leads to 



D (t) - J 2T{t) 1 1 fC7l 

V nm dp (e + 1) 2 (1 + |a 2 ) " ^ 



This is the the expression of the diffusion coefficient D e (t) in eq.(f45"D by the method of the 
Sonine expansion. 

However, the diffusion coefficient D e (t) by the method of the Sonine expansion does not 
work in the Lorentz gas system as follows. From the definition of the diffusion coefficient: 

(r 2 (t)) e = 4 f D e (t)M, (C8) 
Jo 

substitution of eq. flC7|) into this leads to the mean square displacement of the tracer by the 
method of the Sonine expansion as 

<r2(()) - = ^ (TTW^W|^ ' n[1 + ^ po(1 " e2)(1 + l" 2)tl (C9) 

Thus, we have derived the logarithmic time dependence of the mean square displacement in 
terms of a standard method of gas kinetics. 

Figure^ shows (r 2 (t)) e x (1 + e) 2 (l - e 2 )(l + |a 2 ) 2 versus (1 - e 2 )(l + §a 2 )t plots with 
a 2 in eq. (pi8|) for e = 0.6,0.7 and 0.9. Because a 2 for e = 0.4 does not exist, we can not 
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plot the results of our simulation for e = 0.4. It should be mentioned that the results of 



broken line deviates from the results of our simulation. Thus, we suggest that the diffusion 
coefficient D e (t) by the method of the Sonine expansion is not a good approximation in the 
Lorentz gas system. 

APPENDIX D: THE DERIVATION OF DIFFUSION COEFFICIENT D E (V,T) 

From the fact the velocity of the tracer v is distributed, the diffusion coefficient of the 
tracer is expected to depend on the velocity of the tracer v. We derive the expression of 
the diffusion coefficient D e (V,t) in eq. (|48|) which depends on the magnitude of the scaled 
hopping velocity V as follows. 

We derive the diffusion equation which depends on the velocity of the tracer by almost 
similar means which we have adopted to obtain eq. ([4*5]) , i.e. by the means of a Chapman- 
Enskog expansion in the density gradient of the tracer. [[H],|ll[] 

Integrating eq.(|23|) by v , we obtain 



our simulation shown as the symbols are not on a same curve. Equation ( p9| ) shown as the 




(Dl) 



where the density of the tracer n(r, v, t) is defined as 




(D2) 



and the flux of the tracer J(r, v, t) is defined as 




(D3) 



Note that v 



vv with v 



v|. Because eq.( pi| ) holds for each speed of the tracer v, we 



obtain the equation of continuity which depends on v as 



dn(r, v, t) 



+ V-3(r,v,t) = 0. 



(D4) 



dt 
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Here we adopt the Chapman-Enskog method in which the position-velocity distribution 
of the tracer p(r, v, t) is written as a series expansion in a formal parameter e as in eq. . 
Equation fl30|) is now replaced by 

^^f^SfVr.M), (D5) 

01 3=0 

with 

di j) n{r,v,t) = -V (0) ■ J (j) (r,t;,t), (D6) 

and 

3 u) (r,v,t) = J OTp (,) (r. vv,t)6(\v\ - l)dv. (D7) 

Taking into account eq. (^) , the zeroth order of eq. (|24| ) becomes eq. p^ ) . Here eq. (|35| ) is 
replaced by 

P (0) (r,v,t) = !*^/(v,t), (D8) 

where 

m{v,t) = f f(vv,t)5{\v\ - l)dv. (D9) 



From eqs. flD7|) and (|D8|) 



jdO(r iVtt ) =v ^S. [vf(vv,t)6(\v\-l)dv = 0, (D10) 

i.e. the of the tracer vanishes to zeroth order, which leads to d^n(r, v , t) = 0. Addi- 
tionally, from eq.([D^), dp^(r, v, t)/dn(r, v , t) = f(v,t)/m(v,t). 
To the first order in the density gradient e, eq.(^) is replaced by 

{dt n[r ^ t)} dn(r,v,t) om~ 

+ (v ■ V% «, t))^£lM = ApW( r) v , t). (Dll) 
on{r, v, t) 

From the results of the zeroth order and eq.@ with I = l/2dpo, i.e. 
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^ = - V -fh(e)T( t ), (D12) 

eq. (pil| ) becomes 

^h{e)T{t)^ + A)pV>{r,v,t) = (W%M))|^. (D13) 
Similar to eg. (|38|) , the solution (r,v,t) is proportional to V^n(r,v,t) : 



r, v, t) = B(v, T(t)) ■ V (0) n(r, v, t), (D14) 



where the vector B(v, T(t)) also depends on t through T(t). Substitution of this into 



eq. (pi3|) leads to 



(^Me)r (i )^ + A)B(v,T( t ))^^|. (D15) 
Here the contribution to the particle flux of first order in the density gradient is given 

by 

jW(r, v, t) = -D e (v, t) V (0) n(r, v, t). (D16) 

with 

D e (v, t) = - V -Jw- B(w, T(t))6(\v\ - l)dv. (D17) 

Note that we have used symmetry considerations in two-dimension. D e (v,t) is the diffusion 
coefficient which depends on the speed of the tracer v. 

Using this diffusion coefficient D e (v,t) and eq.( P9|) , eq. (|D15|) becomes 



l{>it) He)Ttt)-?- - v D )D e {v,t) = ~. (D18) 



v I w v J dT(t) ' v ' ' 2 

with 



UD ~ fvv-B{vv, T{t))6{\v\ -l)dv ' (D19) 
Dimensional analysis requires D e (v,t) ~ v^ 1 ^) ~ T~*(t). Therefore, 
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dP e (v,t) = D e (v,t) 
dT(t) 2T(t) ' 1 ' 

Substitution of this into eq. (pi8| ) becomes 

v 2 ,h(e)v (t) . i v 2 l _ . 

^M = -(^l +VD ) ->~ Wy (D21) 

where we assume vn ~ vo{t)/l from eg. ( pi9| ). Assuming the hopping velocity 1/r corre- 
sponds to the velocity of the tracer v, we obtain the diffusion coefficient eq. (pE8|) which 
depends on the magnitude of the scaled hopping velocity V as 

D e (V,t) = K^- (D22) 
a + pt 

where k is an unimportant dimensionless constant. 
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FIGURES 




Scatter 

FIG. 1. Unit vector a pointing from the center of a scatter to the center of a tracer. The 
precollision velocity of the tracer and the postcollision velocity of it are expressed as v* and v , 
respectively. 
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FIG. 2. Thermal velocity as a function of time for several values of the restitution coefficient e. 
Note that J(t, e) = (v^it) - w ( 7 1 (0))/(I - e). The straight line 0.058 x tUdr 1 shown as the broken 
line is the best fitting function of the results of our simulation. 
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e 

FIG. 3. M2/2 (short-dashed line) and h(e) (solid line) as a function of the restitution coefficient 
e. M2/2 is obtained from eq, (|B27D , while h(e) is obtained from eq.(||) which has been estimated 
from the data of our simulation. 
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FIG. 4. Mean square displacement of the tracer as a function of time for several values of the 
restitution coefficient e. Note that K(t, e) = (r 2 (t)) e x ln(l + //3). The symbols are results from our 
simulation. The broken line is obtained from eq.(^) with I ~ 13. 5d and ft ~ 0.058 x (1 — e)^" 1 . 
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FIG. 5. Scaled position distribution of the tracer G(R) = \J (r' 2 (t)) e x g e (r,t) as a function of 
the scaled distance R = vj^j (r' 2 (t)) e from the center of the system at which the tracer started. 
The symbols are the data of our simulation which are plotted for several values of the restitution 
coefficient e = 0.4,0.6,0.7 and 0.9 at four different times t = SOOOdU' 1 ^OOdJJ- 1 ^OOdlJ- 1 and 
1000dU~ 1 , respectively. The broken line and the dash-dotted line are the approximate solutions 
of G(R). The simulation data in the tail region correspond to eq.(|58|) with r] = 0.25 and 7 ~ 7.39, 
while the simulation data in the peak region correspond to eq.(|59|) with 77 = 0.75 and V = l/v2- 
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FIG. 6. Scaled velocity distribution of the tracer F(c) as a function of c x . The symbols 
are data of our simulation. The data are plotted for several values of the restitution coefficient 
e = 0.4,0.6,0.7 and 0.9 at four different times t = SOOOdC/- 1 , SOOdlJ- 1 , bOOdU' 1 and lOOOdf/" 1 , 
respectively. The broken line is the Maxwell-Boltzmann distribution 7r _1 e _c2 . 
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FIG. 7. Schematic description of a collision. 
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FIG. 9. Mean square displacement of the tracer by the method of the Sonine expan- 
sion as a function of time for several values of the restitution coefficient e. Note that 
W(t,e) = (r 2 (t)) e x (1 + e) 2 (l - e 2 )(l + |a 2 ) 2 with a 2 in eq. dBlSl ). The symbols are results 
from our simulation. The broken line is obtained from eg . (p9|) with a 2 in eq.( B18| ). 
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